Topological defects and bulk melting of hexagonal ice 
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We use classical molecular dynamics combined with the recently developed metadynamics method 
[A. Laio and M. Parrinello, Procs. Natl. Acad. Sci. USA 99, 20 (2002)] to study the process of bulk 
melting in hexagonal ice. Our simulations show that bulk melting is mediated by the formation of 
topological defects which preserve the coordination of the tetrahedral network. Such defects cluster 
to form a defective region involving about 50 molecules with a surprisingly long life-time. The 
subsequent formation of coordination defects triggers the transition to the liquid state. 
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The melting of ice is a process of obvious universal 
relevance. Under ordinary circumstances melting is nu- 
cleated at surfaces, but situations could be envisaged in 
which the effect of the surfaces is inhibited and melt- 
ing becomes a bulk effect PJ. Studying melting under 
these circumstances can give very precious information 
on the nature of the potential energy surface (PES) and 
on the transition to the disordered state. One could in 
fact imagine that melting takes place either as a sudden 
collapse of the lattice or through the creation and suc- 
cessive catastrophic multiplication of defects . 

In principle molecular dynamics (MD) is an ideal tool 
for studying bulk melting, since the imposed periodic 
boundary conditions eliminate surface effects. Unfor- 
tunately the time scale over which melting occurs is 
too long for present-day computational resources. This 
means that melting can be observed only by super- 
heating the system to the point of inducing a sudden 
lattice instability. Very recently we have developed the 
metadynamics method, which allows long time scale phe- 
nomena to be studied 0, 0, Using metadynamics we 
are able to study the melting of ice near the melting tem- 
perature of the adopted empirical model (^270 K) and 
we find that bulk melting is mediated by the formation 
of topological defects, named "5+7" and characterized 
in a recent work However, in contrast to the picture 
of a sudden multiplication of these defects we find that 
before melting ice goes through a metastable state where 
a defective region of about 50 molecules is surrounded by 
an otherwise perfect lattice. 

In order to overcome the time-scale problem, we ex- 
ploit the extended Lagrangian implementation of meta- 
dynamics 0. The method is based on the construction 
of a coarse-grained non-Markovian dynamics in the space 
defined by a few collective coordinates. The dynamics 
is biased by a history-dependent potential term that, in 
time, fills the free energy minima, allowing an efficient 
exploration and an accurate determination of the free 
energy surface (FES). As in ref. [4|], we introduce an ex- 
tended Hamiltonian that couples the collective variables 



s a (r) to a set of additional dynamic variables s a : 
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where r are the microscopic coordinates of the system 
and Hq the unperturbed Hamiltonian. The masses M a 
and the coupling constants k a should be chosen so as to 
achieve an efficient coupling between the microscopic sys- 
tem and the collective variables, which is obtained when 
the frequencies determined by M a and k a are of the same 
order of magnitude as the characteristic frequencies of the 
microscopic system. The history-dependent potential is 
defined as: 



V(s a ,t) 
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where the time interval At between the placement of 
two successive Gaussians, the Gaussian width Ss and 
the Gaussian height w are free parameters that affect 
the efficiency and the accuracy of the algorithm [jj. 
This method has already been applied successfully to the 
study of rare events occurring in biological systems Q, 
chemical reactions || and phase transitions 0|. 

Two different sets of collective variables were used in 
the study of the melting transition. Common to both 
sets is the use of the potential energy, which is a relevant 
order parameter for the melting as it undergoes a sizable 
change when the phase transition takes place. Moreover 
it is suitable for use as a collective variable in the meta- 
dynamics scheme, as it is an explicit function of the mi- 
croscopic configuration of the system [5j. The potential 
energy was supplemented by coordinates measuring the 
intermediate order. In one case we exploited the Stein- 
hardt |lf| order parameter Qg, which has already been 
used successfully to simulate the nucleation of ice I/, in 
liquid water |ll|. and in the other the number of 5- and 
6-membered rings. In either case we find similar results. 
Here we discuss only results obtained on the basis of the 
second choice. 
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In order to compute the rin g sta tistics we use the short- 
est path ring definition [13, Il3j. This is obtained by 
considering a H 2 molecule and two of its nearest neigh- 
bors (n.n.) and finding the shortest path passing through 
these three molecules. This criterion fails when counting 
large primitive rings but these are of little impor- 
tance for our study. For use in the metadynamics it is 
necessary to turn this definition into a continuous dif- 
ferentiable function of the atomic coordinates Mn({^i}) 
which gives the number of n-membered rings. We de- 
fine the hydrogen bond between water molecules i and j 
by a function of the atomic coordinates, which is 1 
when the two molecules are bonded and otherwise tends 
smoothly to zero For each H2O molecule we con- 

sider the triplets T formed by itself and the pairs of its 
neighbors and we compute the products F T = Y[ t j Cl fij 
of the hydrogen bond functions in the closed paths / con- 
taining T. The shortest path ring containing T is selected 
by a function Gt defined as: 



G T = X/ln e XFT ' /L ' 

IDT 
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where Lj is the topological length of the path I. In the 
limit of large A, Gt tends to the length of the shortest 
path ring. The total number of n-membered rings is then 
given by a sum of rational functions: 



M l ({r 4 })=^- 
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The parameters a and A are chosen so as to achieve the 
correct ring statistics while at the same time making the 
function A/^({rj}) smooth enough to avoid problems in 
the integration of the equations of motion. In our calcu- 
lation we have chosen a = 0.3 and A = 150. 

The initial proton-disordered configurations of ice were 
generated using the Montecarlo procedure described in 
ref. |l6| , which allows supercells with zero dipole moment 
to be produced. We employ the non-polarizable TIP5P 
force field 01 an d treat the long-range electrostatics ex- 
actly by the particle mesh Ewald algorithm. Although 
this model fails in reproducing the thermodynamics of 
water at high temperature (>500 K) [l8|. the hydrogen 
bond dynamics [15j, the self-diffusion coefficient and the 
density anomaly are well reproduced in the range of tem- 
perature of interest to this study (-270 K) [11 EG]]]. 
The melting of ice lh was simulated in the constant pres- 
sure (NPT) ensemble in models consisting of 576 and 360 
water molecules. No size effects on the results were ob- 
served while reducing the size of the system. The results 
reported below refer to simulations of 576 molecules sys- 
tem. In these metadynamics runs the time-dependent 
potential (eqGJ is made up of Gaussians 0.5 kcal/mol 
high, placed every At = I ps. The widths 5s of the 
Gaussians in the space of the collective variables are 10 
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FIG. 1: Evolution of the collective variables during a meta- 
dynamics trajectory for the 576 molecules model. Note that 
the metadynamics time does not correspond to real time. 



and 20 for five and six-membered rings and 24 kcal/mol 
for the potential energy. This choice of parameters leads 
to an accuracy of — 1.5 kcal/mol [20|. 

The evolution of the collective variables during the 
melting is reported in FigQ] Before the phase transi- 
tion takes place, a few events of formation and recom- 
bination of 5-membered rings occur. These events cor- 
respond to the formation of "5+7" topological defects, 
which rccombine after a short time. This defect was re- 
cently discovered by Grishina and Buch j(| and its five 
different conformations arising from the possible proton 
arrangements analyzed. 

In a separate metadynamics run we computed the free 
energy of the defect using as collective variables the ori- 
entation of the O. . .0 bond with respect to two cartesian 
axes , x and z. Simulations were performed in the NVT 
ensemble at two different temperatures, namely 120 and 
270 K. The time-dependent potential is made up of the 
sum of Gaussians 0.167 kcal/mol high, placed every 1 ps, 
with 5s — 0.03 in both the dimensions of the collective 
variables space. 

The FES obtained by metadynamics for the topologi- 
cal defect formation is represented in Fig. |2 The FES 
displays a central symmetry, due to the fact that inter- 
changing the two water molecules used to define the space 
of the collective variable has no effect on the energy of the 
system. Besides the two deep basins (labeled A in Fig. 
I2J) that correspond to the ideal crystal configuration, two 
inequivalent metastable structures have been identified 
(B and C in Fig. They both correspond to 5+7 

defects, but the formation of five and seven-membered 
rings is achieved by different rotations of the pair of H2O 
molecules, leading to structures of varying stability, in 
agreement with the analysis of ref. [f| . The free energy 
of the most stable defect structure (B) is 6.9 kcal/mol 
higher than that of the ideal crystal, which we assume 
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FIG. 2: FES of ice 1^ with respect to the orientation of the 
bond, which is being rotated to form the "5+7" defect at 120 
K. 



as the reference zero value. The free energy barriers for 
the defect formation and recombination are estimated to 
be 8.9 and 2.0 kcal/mol, respectively, which, assuming 
a characteristic frequency for the reaction coordinates of 
~ 5 THz and using the classical transition state theory, 
gives an estimated lifetime of the defect of ~ 0.5 ns at 
120 K. The second defect structure (C) explored during 
the metadynamics run is less stable (F=8.4 kcal/mol) 
and can either recombine or transform into structure B 
through a barrier of 0.9 kcal/mol. Raising the temper- 
ature to 270 K, the height of the free energy minimum 
corresponding to structure B and its formation and re- 
combination barriers remain unchanged, while the recom- 
bination barrier for structure C goes to zero. These free 
energy calculations extend the results of 0] to finite tem- 
perature and confirm the relevance of these topological 
defects also at the melting temperature. 

We found that the "5+7" defects play an even greater 
role, since they are responsible for the shallow minimum 
in the FES indicated as pre-melt in Fig. [3] From metady- 
namics the transition barrier to this local minimum can 
be estimated at roughly 12 kcal/mol. We analyzed the 
nature of this local minimum by performing an inherent 
structure analysis 21] of the metadynamics trajectory, 
quenching to zero K in 50 ps frames of the MD trajectory 
taken every 2 ps. As shown in Fig. the inherent struc- 
tures display a relevant quantity of 5-membered rings and 
a smaller number of 4-membered rings and coordination 
defects. In fact these structures correspond to a con- 
densation of topological defects involving about 50 H2O 
molecules in an otherwise perfect ice Ih lattice. One typ- 
ical defect cluster thus obtained is shown in Fig. The 
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FIG. 3: Two-dimensional projection of the FES into the 
space of the two collective variables energy and 5-membered 
rings. The metadynamics run has been interrupted before the 
basin corresponding to the liquid state was explored. 



energy of the particles, either belonging to smaller rings 
or under-coordinated, ranges from 0.60 to 0.95 kcal/mol 
relative to the energy of ice I/j. 

These defective structures were then embedded in a 
crystalline ice 1^ supercell containing 4608 molecules and 
brought to 270 K, in order to observe their stability. Sev- 
eral MD simulations were performed in the NPT ensem- 
ble with different random initial velocities. Remarkably, 
the average lifetime of the cluster of topological defects 
is 0.4+0.1 ns. This relatively stable accumulation of de- 
fects in a restricted region of the crystal is a nucleus 
for further disorder and melting as the number of co- 
ordination defects and smaller rings grows. In Fig. 0] 
two distinct regimes in the prc-mclting inherent struc- 
tures can be observed. When the system is dragged out 
of the free-energy basin corresponding to the cluster of 
topological defects, a relevant number of small rings and 
coordination defects appears and the energy suddenly in- 
creases. Roughly speaking this signals the watershed be- 
tween configurations belonging to the basin of attraction 
of ice lh and to that of the liquid. 

Since the topological defects cannot migrate, the mo- 
bility of the defective droplet is related only to defect 
formation and recombination at the interface with the 
crystal, and no relevant motion of its center of mass was 
observed. We computed the momentum of inertia of the 
defective region and found that the inertia tensor has two 
eigenvalues T and I2 of similar size and a smaller one I3, 
with an asphericity ratio [/l+j^/^+ja — 0-4- This indi- 
cates that the cluster of defects has a somewhat elongated 
shape. The analysis of the eigenvectors shows that the 
defective region is roughly aligned along the (331) crys- 
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FIG. 4: Defect statistics and energy of the inherent structures 
visited by the system in the pre-melt region. The zero of this 
graph corresponds to the beginning of the melting transition 
in a metadynamics run of a model with 576 water molecules. 




FIG. 5: A cluster of topological defects obtained during the 
inherent structure analysis of the pre-melting phase. H2O 
molecules forming either 4 or 5-membered rings are repre- 
sented by colored sticks. Grey lines represent the ice 1^ struc- 
ture embedding the cluster of defects. 



tallographic direction. Another feature of the FES that 
is worth commenting upon is the shallow basin that ap- 
pears just before the larger liquid basin and corresponds 
to a liquid-solid interface. 

In summary, our simulations reveal that topological 
defects where 5 and 7-membered rings are formed play 
a crucial role in the bulk melting of ice. It is worthy of 
note that in their landmark simulation of ice nucleation 
Matsumoto et al. [22] found that clusters with a similar 
structure to ours provide the nuclei for crystallization. 
The elongated shape of these extended defects and their 
preferential orientation along a given crystallographic di- 



rection might facilitate their experimental detection. The 
recent report that the radiation-induced amorphization 
process in silicon proceeds via the formation and aggre- 
gation of "5+7" defects leads us to believe that this 
is a general feature of disordering processes in tetrahcdral 
networks. 

We thank A. Laio for many useful suggestions, includ- 
ing the use of the potential energy as an order parameter, 
and for reading the manuscript. 
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